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ABSTRACT 

We present a numerical method for studying the normal modes of accretion flows around black holes. In this 
first paper, we focus on two-dimensional, viscous, hydrodynamic disks, for which the linear modes have been 
calculated analytically in previous investigations. We use pseudo-spectral methods and low storage Runge- 
Kutta methods to solve the continuity equation, the Navier-Stokes equation, and the energy equation. We 
devise a number of test problems to verify the implementation. These tests demonstrate the ability of spectral 
methods to handle accurately advection problems and to reproduce correctly the stability criteria for differ- 
entially rotating hydrodynamic flows. They also show that our implementation is able to handle sound wave 
correctly with non-reflective boundary conditions, to recover the standard solution for a viscous spreading ring, 
and produce correctly the Shakura-Sunyaev steady disk solution. Finally, we have applied our algorithm to 
the problem of a non-axisymmetric viscous spreading ring and verify that such configuration is unstable to 
non-antisymmetric perturbations. 
Subject headings: accretion disks, black hole physics, hydrodynamics 



1. INTRODUCTION 

Observations of accretion flows around black holes in 
galactic systems and in the centers of galaxies show strong 
evidence for the presence of long-lived, global modes at 
characteristic frequencies of the black-hole spacetimes (see 
vanderKlis 2000; McClintock, & Remillard 2004, for re- 
views). The orig in of these modes is still a matter of debate 
(see, e.g., [Psaltis 2001, 2004). It is anticipated, however, that 
modeling accurately their origin and physical properties will 
provide the first measurements of the spins of black holes in 
the universe as well as direct tests of strong-field general rel- 
ativity (Psaltis 2004). 

Systematic studies of modes in accretion disks around black 
holes have so far been mostly analytical, performed in the 
linear regime for flows with an artificial viscosity law (see 
Wagoner 1999; Kato 2001, and reference therein; see, how- 
ever, Gammie & Balbus 1994 and Menou 2003 for studies 
of MHD disks). These studies have shown the presence of 
three types of trapped global modes, namely inertial-gravity, 
corrugation, and inertial-pressure modes, with properties that 
appear to agree with observations (see Wagoner et al. 2001 ). 

Despite their success, analytical models have several limi- 
tations. First, they rely on a number of approximations, such 
as the WKB approximation in the radial direction. Second, 
because only linear terms in perturbations are taken into ac- 
count, these studies are unable to reproduce resonances and 
couplings between various modes. However, the detected 
QPOs often have fractional amplitudes as large as 40% of 
the total X-ray fl ux, making the linear approxim ation inac- 
curate (see, e.g., McClintock, & Remillard 2004). Further- 
more, p airs of QPOs are typically observed (see Strohmaver 
2001a b) with frequency ratios equal to ratios of small inte- 
gers, strongly arguing for the presence of resonances between 
modes ( Abramowicz, & Kluzniak 2001). Third, analytical 
and numerical studies of oscillations in MHD disks cast doubt 
on the presence of oscillations at the radial epicyclic frequen- 
cies, as required by all diskoseismic calculations (see, e.g., 
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lHawlev. & KroliiiEooIl 12001 lArmitage. & RevnoldsEool . 
Answering all the above questions necessitates a systematic, 
numerical study of accretion disk oscillations. 

Numerical studies of accretion disk modes have been 
very few and have not addressed in detail the above issues 
(see, e.g., Chen. & Taam 1994, 1995; Milsom. & Taam 1996, 
1997, for simulations of viscous disks; see also Rezzolla 2004 
for simulations of hydrodynamic tori). In fact, no numerical 
simulation of hydrodynamic disks to date has resulted in so- 
lutions that show all three types of accretion-disk modes that 
were discovered analytically. In this series of papers, we per- 
form a systematic study of the oscillatory behavior of hydro- 
dynamic and MHD accretion disks in two and three dimen- 
sions. To this end, we have developed a pseudo-spectral nu- 
merical algorithm, which we present here. Even though spec- 
tral methods have been used extensively in the study of hydro- 
dynamic modes in various planetary atmospheres, including 
the earth's (see, e.g., Mote, & O'Neill 2000), their application 
to the studies of accretion disks has been very limited (Godon 
1997). 

Spectral algorithms are high order numerical methods and 
are perfectly suited to studying modes of hydrodynamic flows 
for a number of reasons. Their major advantage is accu- 
racy and economy in the number of degrees of freedom 
(Bonazzolaet alJ l 19991) . In spectral methods, only ~ n col- 
location points are needed to resolve one wavelength ( Bovd 
2000), whereas in finite difference methods a lot more grid 
points are required to reach the same accuracy. The solution 
is already expressed as a set of orthogonal modes; this prop- 
erty improves the accuracy of studies of mode coupling and 
resonance. Moreover, the dimensionality of the problem can 
be altered very easily, for investigating the effects of the num- 
ber of spatial dimensions on the properties of hydrodynamic 
modes and turbulence. In simulating MHD flows, the high or- 
der of spectral methods allows for the equations to be solved 
in terms of the vector potential and not of the magnetic field it- 
self; this guarantees that the resulting magnetic field is always 
divergence free. Also, super spectral viscosity (Ma 1998a b), 
which is one kind of artificial viscosity, is trivial to implement 



in spectral methods. It resolves shocks and yet preserves the 
high spectral accuracy of the solution. Finally, incorporating 
the self-gravity of the flows is trivial and does not increase sig- 
nificantly the computational time (see, e.g. JChen et alj 2000). 
In this first paper, we discuss our numerical algorithm for 
evolving two-dimensional viscous, hydrodynamic accretion 
disks. In the following section, we present our assumptions 
and equations. In §3, we discuss the details of our numer- 
ical method. In §4, we present a series of tests to ver- 
ify our algorithm. In §5, we apply our algorithm to non- 
axisymmetric instabilities in viscous spreading rings. We con- 
firm the result of Speith. & Kiev ( 2003), which shows that the 
standard solution of a viscous spreading ring (Pi'inglei ri981t 
Frank, King, & Raine 2002) is unstable to non-axisymmetric 
perturbation. In §6, we present our conclusions. 

2. EQUATIONS AND ASSUMPTIONS 

We consider two-dimensional, viscous, compressible flows. 
In this first paper, we neglect self-gravity and the magnetic 
fields of the flows. The hydrodynamic equations, therefore, 
contain the continuity equation 



- + V.( S v)=0, 



the Navier-Stokes equation 



£ — + £(v-V)v = -VP+Vi—Eg, 



and the energy equation 
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We denote by £ the height-integrated density, by v the veloc- 
ity, and by E the thermal energy. In the Navier-Stokes equa- 
tion, P is the height-integrated pressure, r is the viscosity ten- 
sor, and g is the gravitational acceleration. We use $ to denote 
the viscous dissipation rate, q to denote the heat flux vector, 
and F to denote the radiation flux on the r-(j> plane. The last 
term in the heat equation, 2F Z , takes into account the radiation 
losses in the vertical direction. 

We allow for large flexibility in the physical assumptions 
in our algorithm, i.e., we can change the functional form of 
P, g, q, F, and F z easily and apply our algorithm to various 
problems. We also write the viscosity tensor in the general 
form 



Tjj =2(fa + jJ, % )eij+ ( Air + ^b-T/is 



V-v. 



(4) 



Here /x r , /ib, and /i s are the coefficients of radiative, bulk, and 
shear viscosity, respectively. These coefficients can also be 
changed easily to satisfy different physical assumptions. The 
strain-rate tensor e,j is 



1 / dvi dvj 



The viscous dissipation rate is 



(5) 



<i> = 2{ f i r + l i s )(e ij y+U r + fi h --fi s J{V-yY . (6) 

In appendix|X] we write down explicitly the equations used in 
the algorithm. 

The default setup for our calculations is that of a geometri- 
cally thin disk, in which radiation pressure is assumed to be 



negligible. Using the ideal gas law, the thermal energy density 
and pressure become 



£ = £ 



P = E 
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(7) 
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where k B is the Boltzmann constant, T is the central temper- 
ature, fi is the mean molecular weight, and toh is the mass of 
the hydrogen atom. The viscosity coefficients /i r , /ij,, and /i s 
are chosen based on the physical assumptions of the specific 
problems. 

In order to approximate the effects of general relativ- 
ity in the vicinity of compact objects, we use the pseudo- 
Newtonian appro ximation for the gravitational acceleration 
JPaczvnskv. & Wiital 1 9801: iMu khopadhvav 2002) 
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GM 



-2{a/c)y/GMr/c 2 +{a/c) 2 
^/GMr/c 2 (r-r g ) + a/c 



(9) 



Here, c is speed of light, G is the gravitational constant, M 
is the mass of the central object, a is its specific angular mo- 
mentum, r g = 2GM/c 2 is its Schwarzschild radius, and f is 
the unit vector in the radial direction. For a non-rotating star, 
we set a = and the gravitational acceleration reduces to the 
standard pseudo-Newtonian acceleration g — GM/(r— r g ) 2 . 
When a = r g = 0, it recovers Newtonian gr avity. ^_^ 

Regarding radiative cooling, we follow Hubenv ( 1990) and 
Popham, & Naravan (1995) to take into account the energy 
loss in the vertical direction by the standard form 



F 7 = 



3r d 



(10) 



where we use ctsb to denote the Stefan-Boltzmann constant. 
The optical depth r^ is again chosen based on the physical 
problem. We finally set to zero the other two terms, V • q and 
V • F, in the energy equation. 

3. NUMERICAL METHOD 

In our numerical algorithm, we use pseudo-spectral meth- 
ods to evaluate the spatial derivatives in the partial differential 
equations and an explicit Runge-Kutta method to advance the 
solution in time. 

Spectral methods are based on the idea that any function 
can be expanded in a series of orthogonal basis functions. 
When a well-behaved basis is chosen, this series converges 
exponentially in the absence of discontinuities and the nu- 
merical partial derivatives can be easily evaluated. There 
exist different ways to approximate functions by a series 
with a finite number of terms, such as the Galerkin method, 
the ta u method, and the pseudo-spectral method ( see Canuto 
1988; Gottlieb. & Orszagll983HGuoll998tlBovdi20 00: Peyret 
2002). In pseudo-spectral methods, as the one we use here, we 
choose a set of collocation points and then evaluate the expan- 
sion (spectral) coefficients such that the truncated series agree 
exactly with the original functions at the collocation points up 
to the machine accuracy. In the next subsection, we discuss 
our implementation of the pseudo-spectral method. 

3.1. Collocation Methods 

The Fourier basis is the natural choice for expanding func- 
tions with periodic boundary conditions, as is the case along 



the azimuthal direction in our problem. Any physical quantity 
/ = f(r, (j>) can be expanded in this basis as 



/M)= E M r ) e 
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We use a total of M evenly spaced, discrete collocation points 
in the domain [— 7r,7r), 

so that the above Fourier series can be approximated by the 
(truncated) discrete Fourier series 

M/2-1 
m=-M/2 

The derivative of / with respect to is, therefore, 

M/2-1 

= £ imf m {ry 2 ™ 3,M - 

4>j m=-M/2 

In order to evaluate equation J 1 4b . we use a fast Fourier trans- 
form algorithm to find /,„, multiply the Fourier coefficients 
by im, and then take the inverse fast Fourier transform of the 
series. The computational order of the azimuthal derivative is 
therefore &(M\og 2 M). 

When non-periodic boundary conditions are present, as in 
the case along the radial direction, the Fourier basis can no 
longer be used. Popular choices of bases are the Legendre 
and Chebyshev polynomials, which are both defined on the 
domain [—1,1]. Between these two, the Chebyshev polynomi- 
als are usually preferred because of their relation to the co- 
sine function, which makes them easier to compute. Suppose 
the domain of r is [r m i n ,r max ]. In order to map the standard- 
ized radial coordinate r £ [-1,1] to the physical radial coordi- 
nate r £ [r m i n , r max ], we introduce a mapping function r = g(r) 
which is strictly increasing and satisfies both g(-\) = r m j n and 
g(l) = Anax- Let T„(r) be the n-th order Chebyshev polyno- 
mial, i.e., 

T„(r) — cos (narccosr) . (15) 

The Chebyshev-Gauss-Lobatto collocation points are defined 

by 

r k = cos(^Y forO< j<N, (16) 

where the collocation points in the physical radial coordinate 
are given by rj. = g{f] c ). Note that there are totally N +\ 
collocation points including both boundary points of the do- 
main. For any physical quantity f(r, <f>), we approximate its 
Chebyshev-Fourier series by the double sum 

f(rk,4>j)=f{g(h),^j] 

M/2-1 N 

m=-M/2n=0 
M/2-l N 



, COS 



/ TTtlk 



J2iTmj/M 



(17) 
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m=-M/2n=0 \ 

We then use a fast discrete cosine transform to obtain the spec 
tral coefficient /„,„, from /„, (r) = /„, [g(r)]. The radial deriva 
tive can be found using the chain rule 

df 1 df 



dr dgjdr dr 



(18) 



Expressing the derivative with respect to the standardized co- 
ordinate as 






1 N 



dr 



Efe(r) e " 

m=-M/2 n=0 



(19) 



we can employ the following three-term recursive relation to 
obtain /j.'i from /„,„,, 



JN-\,m 

c f (1) 



where cq = 2, and c„ 



fW 

J N,m 



-- o, 

= 2Nf N . m , (20) 

= fnl m + 2(n+l)f n+Lm , 

1, for n > 1. This recursion relation 
allows the use of an &(N) algorithm in calculating the nu- 
merical derivatives in spectral space along the radial direction. 
The transformation between physical and spectral space can 
then be done by a fast cosine transform, the order of which is 
<^(ATog 2 A0. 
The most trivial map r = g(r) is linear, i.e., 

r=—(r + l)- — 



Kr-1). 



(21) 



In this case, 



1 



dr 
dr 



(22) 



dgjdr 

is just a constant. However, the grid spacings Ar^ = r\ — r^_i, 
for 1 < k < N, scale with the number of collocation points Af 
as Ar^ = ff{N~ 2 ), when k is close to 1 or N. This requires the 
stability condition for time-stepping to be 

At = @{N~ A ) (23) 

due to viscous time scale (see SI3.6> . The grid spacings for 
Fourier collocation are A<f>j = <f)j — 4>]-i. They are uniform 
and scale as A0y = &(M~ l ) which gives the stability condi- 
tion 

At = ff{M~ 2 ). (24) 

The stability condition in the radial direction is too expensive 
numerically compared to the one required by the azimuthal 
direction. To overcome this restriction, Kosloff, & Tal-Ezer 
( 1993) proposed the mapping 



^"max 



arcsin (ar) 
arcsin(a) 



+ 1 



' mm 

~2~ 



arcsin ar 



1 



(25) 

arcsin(a) 

which makes the spacing Ar/, — 6[N~ ) around the bound- 
aries. In the above transformation , a is a paramete r. Let e 
be the machine accuracy. Don, & Solomonofn i 1 9971) showed 
that the choice 

""^ ,26, 



sech 



N 



minimizes the errors in the numerical derivatives. We there- 
fore refer to the mapping d25l > with the optimized choice J26I 
as the modified Chebyshev collocation and employ it in our 
algorithm. The numerical derivatives for the modified Cheby- 
shev collocation can be easily calculated by 



df 
dr 



2arcsin(a) 



('"n 



a) a 



1 ( \2 9 f 

1 " {ark) 1 ^r 

or 



(27) 



where (df/dr)\ ?k is obtained by the usual Chebyshev 
method dl9> . The constant scaling in (I22t and the scaling in 
equation ( I27> can be done together with the normalization, so 
neither of them have extra computational cost. 



3.2. Treatment of the Non-Linear Terms 

In low order numerical methods, one uses the conservative 
form of the hydrodynamic equations to ensure conservation of 
mass and momentum (LeVeque 1992). This is not necessarily 
the ideal way for pseudo-spectral methods because t he con- 
servative form is sometime numerically unstable (see Pevret 
2002, pp. 286-289). To illustrate this, consider two arbitrary 
functions u(x) and v(x). The derivative of their product has 
two different but equivalent forms 



C=(«v)', 

D = U V+UV . 



(28) 
(29) 



The pseudo-spectral Fourier coefficients (with M colloca- 
tion points) of C and D are 



C„ = im Y, u p v q + i 

p+q=m 

+ Y, rnu p v c 

p+q=m-M 



D m = im y\ " 

p+q=m 



pVq + l 



£ rnu,pV q 

p+q=m+M 



V {m+M)u p Vq 

p+q=m+M 



(30) 



+ £ (m-M)u p v q 

p+q=m—M 



(31) 



The first term in each equation corresponds to the truncated 
Fourier series. The other terms in the square brackets are 
ther efore due to t he aliasing error of pseudo-spectral methods 
fsee lBovdl 12000. pp. 202-221), which may cause numerical 
instability. One can show that these aliasing terms have oppo- 
site signs. The skew symmetric form (C m +D m )/2, therefore, 
has a much smaller aliasing error. Unfortunately, calculating 
the skew symmetric form increases the number of numerical 
derivatives and lowers the algorithmic performance. 

A similar argument is valid for the Chebyshev colloca- 
tion. However, in that case, the problem becomes more com- 
plicated due to the boundary conditions. Botella. & Pevret 
(2001) carried out numerical experiments and reported that, 
without aliasing removal, the convective form is stable for the 
two different boundary method they tested; the conservative 
form is unstable for one of their methods. We carried out 
some numerical experiments for our algorithm. When the ini- 
tial conditions are smooth, both conservative and convective 
forms are stable. However, when an initial perturbation is in- 
troduced, the conservative form sometimes become unstable. 
In these cases, the spectral coefficients at high wavenumbers 
in the radial direction grow exponentially even with spectral 
filtering (see the next subsection). The convective form, on 
the other hand, is able to reproduce various analytical stabil- 
ity criteria, as shown in |4] 

Because the convective form gives an approximation that 
is of the same order as the conservative form, in the present 
implementations of the numerical algorithm, we use the con- 
vective form for the non-linear terms in Navier-Stokes equa- 
tion to ensure numerical stability. The continuity and energy 
equations, on the other hand, are written in conservative form 
to conserve mass and energy. The exact equations we solve in 
the algorithm can be found in appendixlAl 

3.3. Spectral Filtering 



In our numerical algorithm, we use a combination of spec- 
tral filters in order to resolve the two principle drawbacks of 
spectral methods. First, because of the non-linear character 
of the Navier-Stokes equations described above, we need to 
filter out the high-frequency modes in each time step to re- 
duce the aliasing error. The spectral filters ensure long-time 
stability of the algorithm. Second, when shocks are present 
in the solutions, the spectral coefficients do not converge ex- 
ponentially and oscillations appear around the discontinuities 
(Gibbs phenomenon). Introducing an additional spectral fil- 
ter increases the converge rate of the solutions. Indeed, it has 
been proven that the filtered solution converges to the correct 
entry solution and hence gives the correct sh ock properties. 
(For more information see|Don||T 994; Don, & Ouillen 1995; 
iGottlieb. & Shull997tlMall998albHGuo et al.l200lULil200ll 
and reference therein.) 

Suppose we want to apply a filter in the (^-direction. We first 
find the Fourier coefficients f m (r) according to equation (lilt . 
We then denote the filtered sum by 

f v (r,<t, J )= £ <rJ — )f m (r)e 

m=-M/2 V 1V1 / 

where we use the exponential filter 



ilTtmj jM 



(32) 
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exp — | Inej 
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M 



(33) 



The parameter e here is again the machine accuracy. We 
choose j3 ~ M/2 so that the filtered sum approximates the 
original function very well and does not reduce the accuracy 
of the numerical solution. The same filter is used in the radial 
direction, i.e., 



0/3 



exp 



I In el 



(34) 



where j3 is chosen close to N so the double filtered sum is 

M/2-1 N 
m=-M/2n=0 VV 

/ imk\ 
. I 



2m\ 



Xfn,, 



V N 



M ) 



3.4. Initial Perturbations 



Spectral methods have been used successfully to study non- 
randomized perturbation, such as the tidal effects presented 
by Godon ( 19971 11998ft . Later. iGodon & L ivio (2000) de- 
scribed a method to add random perturbations to their ini- 
tial conditions. In their simulation, they generated 100 Gaus- 
sian profiles randomly within the computation domain. This 
method ensure that the spectral coefficients of the perturbed 
function converge exponentially, and avoids changing the di- 
rection of the characteristic at the boundary. 

In our calculations, we require the initial perturbations to 
be smooth, i.e., we require the spectrum of the perturbations 
to converge exponentially. One simple method to obtain such 
a perturbation is to generate random noise in spectral space, 
multiply the noise by an exponential filter, and transform the 
noise back to physical space. 

To formulate this method, we first generate a flat-noise 
spectrum, P njn , which is totally random. The perturbation is 
added to the original function by the following equation 



/M) = [l+7(r,<^^M)]/oM), 



(36) 
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FIG. 1. — Two 1 -dimension perturbation profiles generated by the method describe in 33.41 All lines in the plots have 70 = 0.001. (Left) The function is assumed 
periodic and approximated by a Fourier collocation. (Right) A Chebyshev collocation is used to demonstrate the behavior of our method near the boundaries. 



where f(r,c/>) is the perturbed function, j(r,(f>) is the size of 
the perturbation, P^P^P (r, cj>) is the filtered sum 

M/2-1 N f v 



m=-M/2»= 



X Zjm ^^ CO s 



/? 



/™£\ Uirmj/M 



(37) 



as defined in the pervious subsection, and fo(r,4>) is the un- 
perturbed function. The order of filter, (3, controls the spacial 
scale of the perturbation. Because fo(r,<j>) is periodic in <fi, 
we can simply use a ^-independent function j(r,(j)) = j(r). 
However, in the r-direction we do not want the perturbation 
to extend the boundaries. We, therefore, set 



l{r) = 70 



arcsin(a) 



l-(ar) 2 , 



(38) 



which is pre-computed by equation (I27> . when we use the 
modified Chebyshev collocation method. 

In Figure [Q we illustrate this method of introducing per- 
turbations to an one-dimensional constant function /o = 1. 
The left plot uses a Fourier basis with 256 collocation points. 
Where the right plot uses a Chebyshev basis with 257 col- 
location points. All the lines have 70 = 0.001 but different 
orders (3. Note that the extra rescaling for Chebyshev collo- 
cation lowers the amplitude of the perturbations around the 
boundaries. 

3.5. Boundary Conditions 

Solving partial differential equations in a finite domain usu- 
ally requires Dirichlet or/and Neumann boundary conditions. 
When using finite difference methods, one can simply change 
the values of physical quantities at the boundary points to 
achieve Dirichlet boundary conditions. Alternatively, one can 
adjust the "ghost points" to achieve Neumann boundary con- 
ditions. Unfortunately, these methods do not generally work 
for spectral methods. Indeed, changing the boundary points 
(or ghost points) effectively introduces step functions. When 
discontinuities are present, the spectral coefficients do not 
converge exponentially and the spectral methods fail to pro- 
duce stable solutions. 

Perturbations around the boundaries can also change the 
characteristic directions of the flow. A naive Dirichlet or Neu- 
mann boundary treatment fails to take care of these oscilla- 
tions between inflowing and outgoing characteristics. When 



the inflowing boundary conditions are not given, or when the 
outgoing boundary conditions are not consistent with the so- 
lutions in the computation domain, the system of differential 
equations is ill-posed and the algorithm diverges. Although 
we generate the initial perturbations in a way that they are 
small around the boundaries (see the previous subsection), 
they can still propagate and reach the boundaries as time 
evolves. 

In our numerical algorithm, we introduce a new boundary 
treatment in the radial direction, which is a spatial filter. This 
method forces each dynamic variable to approach smoothly 
to its boundary value so the characteristics directions at the 
boundaries are well-defined. This ensures that any instabili- 
ties of the solution are not due to boundary conditions. For 
example, consider a one-dimensional problem 



df(r,t) 
dt 



■F\f{r,t)\ 



(39) 



for any functional form F. If we want to apply a Dirichlet 
boundary condition at the outer boundary, we choose some 
smooth monotonic function h(r) such that h(r) — > 1 for inte- 
rior points of the computational domain and h(r) — > at the 
boundary. At each time step, we impose the boundary condi- 
tion by the spatial filter 



fl-^(fl-fo)h(r k )+f , 



(40) 



where fi denotes the numerical approximation of our func- 
tion at time i and radial collocation point k. This filter makes 
the function / approach its boundary value /o and produces 
a numerical boundary layer. Rearranging, the above step is 
equivalent to setting 



fl-*fl-[l-h(r k )}(fl-f ). 



(41) 



Hence, it is equivalent to adding an extra source/sink term / s 
in the original equation 



df 

dt 



F[fV 



l-h(r) 
At 



(f-h) = F [/]+/;. 



(42) 



where At is the time step (see 33.6J . Omitting F in the 
above equation, it becomes df /d ~ _ (/ _ /o)- The value of 
[1 —h(r)]/ At controls the converge rate of / to /q, which is au- 
tomatically proportional to the velocity, because 1/Af ~ |v|. 



We use an exponential filter that is similar to the spectral fil- 
ter discussed in §3.3. To impose an outer boundary condition, 
we need the monotonically decreasing filter 



h{r) = ap 



max ' mm 



= exp 



I In el 



max ' mm 



(43) 

On the other hand, for the inner boundary, we use the mono- 
tonically increasing filter 



h{r) = ap 



max ' mm 



exp 



I In el 



max ' mm 



(44) 

We can change the thickness of the numerical boundary layer 
by changing the order of the filter, f3. 

This boundary treatment turns out to be very convenient 
in modeling non-reflective boundary conditions. When we 
choose the parameter (3 to be equal to the number of colloca- 
tion points, the boundary layer is thick enough to discard any 
outgoing waves but thin enough so that the interior solution is 
not affected significantly. 

For applying Neumann boundary conditions, the addition 
of a ghost zone is not practical in spectral methods. This is 
because Chebyshev polynomials approxi mate th e derivatives 
based on all the interior points. Godqn| (119971) presented a 
method which involves solving for the boundary values of 
physical quantities, so that their derivatives agree with the 
boundary conditions. This is effectively a Dirichlet bound- 
ary treatment with the boundary values depending on all the 
interior values at each time step. 

To avoid numerical instabilities due to the Neumann bound- 
ary conditions, (<9//9r)|t,oundary = /n, we choose another 
method that involves applying a spatial filter to each variable 
after taking the radial derivatives. This is equivalent to replac- 
ing the radial derivatives (df/dr)\,- k by 



df 

Or 



df 

Or 



1-0)9 



n-r B 



max ' mm 



-fo 



for the inner boundary, and 



(45) 



df 

Or 



df 

Or 



1-0)3 



- n< 



max ' mm 



■fo 



(46) 

for the outer boundary. This ensures that the first derivatives 
around the boundaries are smooth. 

3.6. Time Stepping 

We integrate the hydrodynamics equations with low- 
storage, explicit Runge-Kutta methods. Let 

U=(H^H Vr ,H^,H E )' (47) 

be th e right hand sides of hydrodynamic equations JA1I - 
(lA4l and let 

u = (E,v„v*,£)' (48) 

be the dynamic variable; the superscript t here stands for 
transpose. The low-storage explicit Runge-Kutta methods 
have the general from 

u = u(r„) 

Q, = AjQm+A*H(«m) 
Uj = u,_i+B,Q,-, i= l,...,s ' 
u(f„+i) = u s 



(49) 



where At is the time step and A\ = 0. All other A,'s and 
Bi's are coefficients that characterize the scheme and u, are 
the intermediate stages of u. The method to obtain these co- 
efficients can be found in Pevret (2002). Depending on the 
complexity of the problem, we choose between second order 
or third order Runge-Kutta, and third order or fourth order 
Carpenter-Kennedy Runge-Kutta methods. The different al- 
gorithms are compared during the verification tests (see 34.11 
and fig.|3]below). The third order Runge-Kutta scheme, pro- 
posed by Williamson (See Pevret 2002, p. 146), turns out to be 
the most popular in spectral methods because of its efficiency 
and accuracy. For this reason, we present here this scheme 
explicitly: 



u« 


= u(f„) 


Q: 


= AfH(uo) 


Ul 


= U0+5Q1 


Q2 


= -§Qi+A*H(u,) 


U 2 


= "1 + t!Q2 


Q3 


= -|§Q 2 + AfH(u 2 ) 


U 3 


= "2+^Q 3 


(t„+l) 


= U 3 . 



(50) 



Recalling the definition of Ar^ and A<pj in 33.11 we de- 
fine the Courant-Friedrich-Lewy (CFL) time step in different 
directions independently, 



M 



CFL,r = min 



C s + V, 



AfcFL,0 = mm , (51) 



where c s = y/P/Tl is the sound speed. Let Al = 
min(Art,rtA0j) be the minimum grid separation. We also 
define the viscously restricted time step as 



At,, 



AV 



(52) 



where u s — /x s /£ is the kinematic viscosity coefficient. Us- 
ing these definitions, the maximum allowed time step At is 
chosen adaptively by 



At = (5min(Af C FL. r , Af C FL,0, Af„ 



(53) 



Here 6 is a constant that depends on the physical problem. For 
example, for a free falling dust ring, the flow is very stable 
and 5 could be chosen as large as 29. For most of accretion 
problems with vy <C Va,, S is of order unity. 

4. CODE VERIFICATION 

We have verified our numerical algorithm by comparing our 
numerical results to a number of test problems with analytical 
solutions. These problems are designed to test the implemen- 
tation of the different terms in the equations, such as the ones 
that describe gravity, pressure, and viscosity. In this section, 
we describe these tests in details 2 . 

4.1. Free Fall of a Dust Ring 

As a first test, we study the free-fall of an axisymmetric dis- 
tribution of matter with no angular momentum, pressure, or 
viscosity. This problem tests the implementation of gravity, 
the radial boundary conditions, and the conservation of mass 

2 The simulation corresponding to these test can be find in the webpage 

http : //www. physic 5 . arizona . edu/ ~chan/research/ 
|astro-ph/ 040 6073V . 
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in the radial direction. Also, we use this problem to evalu- 
ate the efficiency and stability of the different time-stepping 
schemes. We set v^ = 0, P = 0, i/ s = 0, and consider only 
Newtonian gravity with g = (GM/r 2 )r. 

In order to verify our numerical result, we also solve the 
problem analytically. Using dimensionless quantities so that 
GM = 1, we write the conservation of energy for a fluid ele- 



ment initially at tq as 



drY 

dt ) 



2 



(54) 



Integrating over time, we obtain the trajectory r = r(t;ro) of 
each fluid element in the implicit form 

^ , ! • (, F 

— TTT-f = - sin 2 arccos , / — 



+ arccos 



(55) 



8 



We now consider conservation of mass, i.e., at any time we 
require that 



2irrY,(r,t)dr = 27rroS]o('"o)^''o- 
Evaluating dro by using equation ( 1551 . i.e., 



dr = 



rodr 



3/2 r 

2 V r y/ro/r- 1 



the analytical solution for the density becomes 

,•2 



E(r,f) = £ fa,)-£ 




-,-1 



(56) 



(57) 



(58) 



where ro can be obtained implicitly using equation d55i . 

Although this is an one dimensional problem, we simulate it 
in the two-dimensional domain [0.2, 1.8] x [— 7r,7r] with 257 x 
64 collocation points. The initial condition is a Gaussian ring 

:exp[-20(r-l) 2 ] 



£o(r,( 



(59) 



with zero initial velocity. At the outer boundary, we use 
the standard Neumann boundary conditions, i.e., we set 
(9S/5r)| roul = and (dv r /dr)\ rout = for all times. More- 
over, because the characteristics at the inner boundary point 
towards the negative r-direction and we have neglected vis- 
cosity and pressure, we do not need to impose any explicit 
boundary conditions there. Because filtering stabilizes the 
flow, we do not apply any filters in this test in order to com- 
pare how different time-stepping methods affect the stability 
of the algorithm. 

The top panels of Figure |2]compare the results of the simu- 
lations to the analytical solutions; the open circles denote the 
numerical results using the third order Runge-Kutta method 
and the solid lines are the analytical solution. The numerical 
results are indistinguishable from the analytical solution. The 
bottom panels of Figure|2]show the error between the numer- 
ical and analytical solution for the maximum stable time step. 

We test four different time-stepping schemes: second order 
Runge-Kutta (RK2), third order Runge-Kutta (RK3), third or- 
der Carpenter-Kennedy (RK3CK), and four order Carpenter- 
Kennedy (RK4CK). Figure |3] shows the result of this test. In 
the left plot, the vertical axis shows the computation time. The 
horizontal axis shows the dimensionless t ime- stepping factor 
5 = Af/min(Af C FL : r,AfcFL.0,A/ l/ ) (see J3~6l The right plot 
shows the numerical errors, max(|S num -E ana |), at t = 1, for 
the same problem. For 6 > 15 in RK2, 6 > 18 in RK3, S > 21 
in RK3CK, and S > 29 in RK4CK, the solutions diverge so 
neither the computation time nor the error are plotted in the 
graphs. This test demonstrates the ability of our algorithm 
to conserve mass to an accuracy better than a few parts in a 
million. It also demonstrates the advantage of the third order 
Runge-Kutta method over other algorithms, when solutions of 
high accuracy are necessary. 

4.2. Rayleigh 's Criterion 

We now study the ability of our algorithm to model cor- 
rectly instabilities in hydrodynamic flows. In the absence of 
viscosity, the necessary and sufficient condition for a rotating 
flow with angular velocity f2(r) to be stable is 



-L(r 2 n) 2 >0 7 
dr 



(60) 



which is known as Rayleigh's criterion ( Chandrasekhar 
1981). Although Rayleigh's criterion usually describes the 



TABLE 1 

Numerical Sound Speeds for the Problem discussed in §4.4. 



r 


Numeric Sound Speed 


Analytical Sound Speed 


Percentage Error 


1.0 


1.011210 


1.000000 


1.12095% 


i.i 


1.062727 


1.048808 


1.32702% 


1.2 


1.107763 


1.095445 


1.12449% 


1.3 


1.153139 


1.140175 


1.13695% 


1.4 


1.193408 


1.183216 


0.86140% 


1.5 


1.235294 


1.224744 


0.86136% 



stability condition in a Couette flow, i.e., for an incompress- 
ible fluid between two rotating cylinders, the criterion is also 
valid for rotating compressible fluids, when the unperturbed 
radial velocity vanishes everywhere in the flow. This can be 
achieved by balancing the centrifugal and centripetal acceler- 
ations. Let 



/ dfl 

K 2 = 2n\2n + r — 

dr 



(61) 



be the square of the epicyclic frequency. Rayleigh's criterion 
is equivalent to requiring n 2 > 0, because d(r 2 Vt) 2 / dr — r 3 n 2 . 
We define a modified Newtonian gravity with a negative 
gravitational index a < such that the gravitational accelera- 
tion is 

*(r) S -^. (62) 

r\ a \ 

Setting GM = 1, the corresponding "Keplerian" velocity v^ is 
given by 

V4> = Vr^. (63) 

Substituting this into Rayleigh's criterion, i.e., equation J60l >. 
we find that a — -3 is the critical stability condition. The 
growth rate of the instability is 



= Im-v/ft? = ImJ (a + 3)r a ~ l 



(64) 



Of course, the same result can be obtained by observing that 
the effective potential 



2r L \Oi+lj 



v+l 



(65) 



has no local minimum when a < -3; here L is the angular 
momentum of a fluid element with a unit mass. 

We use a spatial resolution of 257 x 64 collocation points 
to test the Rayleigh criterion. The initial conditions are with 
a uniform density, So = 1, and the velocity profile (I63> . Ran- 
dom perturbations that are of order 10~ 6 are added to all 
physical quantities and the boundary conditions are set to 

^Kmin = V n>max = "■ 

Because Rayleigh's criterion refers to an one-dimensional 
problem, we discuss here explicitly the evolution of the m = 
mode. This is equivalent to averaging over the azimuthal di- 
rection. Figure@]shows the evolution of power at k = 16, 32, 
64, and 128 (and m = 0) of the density, i.e., it is the plot of 
|£fc m (?)| 2 against time. When the flow is unstable, the pertur- 
bations grow very fast and diverge. As required by Rayleigh's 
criterion, flows with a > -3 are stable to perturbations. 

4.3. Propagation of Wavefronts 

In order to test our implementation of the thermal pressure 
and the numerical spreading of wavefronts, we neglect grav- 
ity, set all velocities equal to zero, and simulate the propaga- 
tion of a sound wave. We use the equation of state P = KT, r 
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FIG. 4. — The evolution of the power at modes with k = 16, 32, 64, and 128 (and m = 0) of the density for the study of Rayleigh's criterion discussed in i|4.2l 




FIG. 5. — The propagation of a sound wave in a uniform static background. The gray-scale snapshots correspond to the fluid density; dark is high and light is 
low density. 
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FIG. 6. — The evolution of a viscous spreading ring. The solid line is the 
analytical solution and the circles are numerical solutions. Note that con- 
centric rings appear at the beginning, but, after about 24 ring rotations, they 
decay away and we recover the analytic solution. 



with K=\ and vary the polytropic index T from 1 to 1 .5. This 
pressure term decouples the energy equation from the Navier- 
Stokes equation and the sound speed is simply vaT5?^. 
The initial density is 



Eo(r,( 



l + 10" b exp -60(l + r /! -2rcos( 



(66) 



We use 257 x 64 collocation points and calculate the solution 
from t — to 2. 

Figure |5] shows six gray-scale snapshots of the density for 
T = 1 . The resulting sound wave propagates outward as a cir- 
cle. Note that, because of our boundary treatment, the density 
wave is not reflected back as it propagates close to the bound- 
aries. To measure the sound speed, we take the density profile 
at (p — and trace the wavefronts. The collocation points that 
correspond to the peak density are then used to calculate, via 
X 2 -fitting, the sound speed. Table^shows the x 2 -fitted sound 
speed of the simulation. Compared to the analytical sound 
speed, the error is of order 1%. 

4.4. Viscous Spreading of a Fluid Ring 

The spreading of an axisymmetric viscous ring with small 
kinematic viscosity v s (lPringldll981t iFrank. Kin g. & Raine 
2002) has become a standard test for numerical methods in- 
cluding both smoothed particle hydrodynamics (SPH) and 
grid-based codes. It has also often been used to guide studies 
of many other astrophysical problems, su ch as the fall-back 
disk after a supernova explosion (see, e.g., van Paradiis et al. 
1995; Eksi.&Alpar 2003). 

In the standard solution, it is assumed that v^ = rilx 3> v r , 
so the viscosity tensor can be approximated by 



= 0, T r< j> — T^r = V S Y> 



dv<f, 
dr 



V£ 

r 



(67) 



Using conservation of mass, conservation of angular momen- 
tum, the fact that v^ is always close to the Keplerian value, 
and assuming that the kinematic viscosity coefficient v s is a 
constant, the analytical solution for the density, in our units, 
becomes 



E(r,0 



1 



YI-kv^IH 



exp 



l + r z 
\2vj 



'1/4 



2r 
\2v~t 



(68) 



where I\ u is a modified Bessel function. 



We perform our s imulations with axisym metric initial con- 
ditions. We follow Speith, & Kiev] d2003l) and use a typical 
value v s — 4.77 x 10~ 5 . The energy equation is decoupled 
from the hydrodynamic equations as in the previous subsec- 
tion. We set r = 1, and choose K such that c s = 10~ 8 . The 
spatial resolution is 257 x 64 and the time step factor 6 = 0.8. 
We start with an initial density given by equation fl68l > at 
t = 27.95, which corresponds to 12^ s f = 0.016. The boundary 
conditions are v |,- max = (rtt K )\ rmx and v | rmin = (r^)| rmin , 
which are the ones assumed in the standard solution. 

Figure [6] compares the analytical solution to the numerical 
solution. The simulation shows that some concentric rings ap- 
pear at the beginning. They come from the relaxation effect 
as our initial velocities approach standard solution. In fact, 
these additional structures have been found in other simula- 
tions (for more information see Monaghan 1992), and decay 
away after about 24 rotations. 

4.5. Shakura-Sunyaev Steady Disk Solution 

Here we consider the thin disk problem and suppose that 
the disk is able to settle to a steady-state structure. By assum- 
ing an a-viscosity law, Shakura, & Sunvaev ( 1973) solved the 
local disk structure analytically. Their work has been highly 
cited for studies in accretion disks. (See Frank, King, & Raine 
2002 and reference therein). 

In our simulation, we solve the hydrodynamic equations Q 
- including the energy equation. The kinematic viscosity 
coefficient is given by 



= ac<;H. 



(69) 



The sound speed is given by cl = dP/dY, — k^T / ' pmi^ and the 
disk scale height is H = rc % jv^. We, also, assume the vertical 
optical depth of the disk to be t& — Ekr, where the Rosseland 
mean opacity kr is given by Kramers' law 3 

E 



kr = 6.07 x 10 



22 



-7/2. 



A-kH 



(70) 



Based on the above assumption and using the fact that u s is 
small, the Shakura-Sunyaev disk solution is given by 



E = 6.15 «- 4 / 5 < /IO m! /4 ^ /4 / 14/5 g cm 



v, = -2.59xl0V/ 5 < /10 m7 1/4 «>- 1/4 - 



K \Q J 



14 / 5 cm S -' 



v = 1.15x I0 & m\ /2 R 



M 



-1/2 
10 C 



T=1.48x 10V 



■'^V 4 * 



"3/4 f 6/5 K 
10 / ^' 



(71) 
(72) 
(73) 
(74) 



where M 16 =M/10 16 gs -1 , mj =M/M & , R l0 = r/(10 10 cm), 
and / = (1- y/R^/r) 1 ' 4 , with R* being the radius of the cen- 
tral object. Considering the flow around a standard neutron 
star with a mass of 1.4M , and a radius of 8.5 x 10 5 cm, 
we setup the simulation with r m ; n = 10 9 cm and r max = 1.2 x 
10 n cm. The initial conditions for density and energy are as 
the Shakura-Sunyaev disk solution with M\s = 1, ni\ = 1.4, 
and a = 0.1. We set the initial radial velocity to zero and 
the initial azimuthal velocity to be Keplerian. Perturbations 
were added to all variables with amplitude 0.1 and order 3. 
The simulation runs with 513 x 64 collocation points up to 
t = 10 7 s, which is approximately 16 viscous time scales. 

In Figure we show our numerical steady solution for 
a = 0.1, atf = 10 7 s; note that we use the magnitude of the 



3 See again Frank. Kine. & Raine i2002), p. 94 for notes on errors in the 
literature. 
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FIG. 7. — (Top) The analytical solution (solid line) compared to the numerical solution (open circles) of the Shakura-Sunyaev steady disk solution in log-log 
scale. (Bottom) The percentage numerical errors, defined as (fnum/fimn _ 1) x 100%. All the units along the vertical axis are cgs units, i.e., g cm" 2 for density, K 
for temperature, and cm s for velocities. 



radial velocity in order to plot it in log-log scale. The numer- 
ical values are different from the Shakura-Sunyaev solution 
by a few percent. We point out explicitly that the pressure 
term in the Navier-Stokes equation is neglected in obtaining 
the analytical solution. Although the pressure term does not 
transport angular momentum, the effect of pressure compared 
to those of viscosity is proportional to I /a. Since a — 0.1 in 
our simulation, the pressure effects can account for the differ- 
ence between the analytical and numerical solution. 



APPLICATION: NON-AXISYMMETRIC INSTABILITIES IN 
VISCOUS SPREADING RINGS 



5. 



As we mentioned in 34.41 the viscous spreading ring has 
become a standard test for numerical simulations. In these 
simulations, additional structures always appear, such as the 
concentric rings in Figure |6] We attributed the presence of 
these concentric rings to relaxation effects of the initial con- 
ditions. Maddison et al. (1996), on the other hand, showed 
that the concentric rings that appear in SPH simulations using 
the artificial viscosity given by Monaahan ( 1992) are numer- 
ical artifacts. Later, Speith, & Kiev (2003) used perturbation 
methods to show that viscous spreading rings are unstable to 
non-axisymmetric perturbations. 

Assuming a compressible, pressureless, and viscous fluid, 
the equations used in Speith. & Kiev ( 2003 ) are the continuity 
equation ([0 and the Navier-Stokes equation (|2j, with zero 
pressure. By setting v % = /j, s /H to be constant and /i r = /j, b = 
0, the first order instability is an one-armed spiral, whereas 



the second order instability is a superposition of an one-armed 
and a two-armed spirals. 
If we define a function 



E(t r)-2-^t 



(75) 



where 51k = \/GM/r i is the Keplerian angular velocity, then 
E(t,r) measures the first order perturbation amplitude. Its dis- 
persion relation is 



( v s 8v s <9£ 



+i r- 



2"s 



v s 9Eq 



3 r£ dr 



(76) 



where k is the radial wavenumber and So is the standard so- 
lution J68> . Assuming llvj <C r, it follows that the condition 
for instability is 

3 dt ° HI) 



k 2 > 



rS dr 



and the growth rate is 






(78) 



In order to test this claim, we carried out a test with the same 
condition as Speith and Kley's RH2D simulation (a radiative 
hydrodynamic code; see Kiev 1989). The only differences are 
in the initial perturbations and a small constant bac kgrou nd 
term in the density, i.e., the setup is the same as in 34.41 but 
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FIG. 8. — Snapshot of simulations of a viscous spreading ring with different resolutions. Note that for different resolution, the direction of the spirals are 
different. 
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FIG. 9. — The evolution of a viscous spreading ring with non-axisymmetric perturbations. At t = 100, there are concentric rings due to the initial relaxation. At 
t = 200 it is clear that a small, one-armed spiral starts to develop in the inner region. The remaining figures, then, show the development of the one-armed spiral. 
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FIG. 10. — The contour plot of density at (j> = against time. Region I 
shows the relaxation effect. Region II shows that the solution has relaxed 
to the standard solution. The Speith-Kley instability develops in region III. 
In region IV, the unstable wavelengths are too short to be resolved in our 
numerical method. 
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FIG. 11. — The evolution of the first 6 modes of E(r = l,f). 



we add initial perturbations by the method described in 33.41 
with order 3 and magnitude 0.001. The additional constant 
background term in the density varies from 0.1 to 0.001 de- 
pending upon the resolution, which keeps the solution stable. 
We perform the simulation with resolution 129 x 64, 257 x 64, 
513x64, and 1025x64. 

Figure [8] shows the density contours for different resolu- 
tions at time t = 120. It is clear that the direction of the spiral 
in low resolution is different then in high resolution. We con- 
firm this result with Speith and Kley (private communication), 
i.e., that the transition radius between leading and trailing spi- 
rals vary as the resolution changes. Figure [9] is the density 
contour of the simulation with 513 x 64 collocation points 
and 0.001 background density, at different times. At t = 100, 
some concentric rings due to relaxation effects are seen. At 
t = 200, a one-armed spiral starts to form at the inner disk. 
The spiral is fully developed by t = 300. 



In order to find out the unstable wave number, we plot in 
Figure^|the contour of density along = against time. Re- 
gion I shows the relaxation effect, where the relaxation waves 
propagate outwards. In region II, the solution has relaxed to 
the standard solution. Later on, the Speith-Kley instability 
starts to develop, which is shown in region III. We can es- 
timate the wavelength of the fully developed spirals, which 
varies from 0.1 to 0.15. For region IV, the unstable wave- 
lengths are too short to be resolved in our numerical method. 

Finally, in Figure ^2 we pl°t the first 6 modes, i.e., m = 
0, 1,. . . ,5, of E(r = l,t). The late time behavior of our solu- 
tion agrees with Figure 8 in Speith, & Kiev (2003). To com- 
pare our numerical result with the analytical growth rate J78t . 
we X 2 "fit the slope of the m = 1 mode in Figure fTTIbetween 
t = 250 and t = 500, which gives Im(er d ) sa 0.0013. That cor- 
responds to unstable wavenumber k sw 9.0. It agrees with the 
width of the spiral ss 0. 1 1 as shown in Figuresl9landllOI 

Physically, the viscous spreading ring problem is an im- 
portant guide to understanding the transport of angular mo- 
mentum in accretion disks. We have independently verified 
the Speith-Kley instability: a cold viscous disk is unstable to 
non-axisymmetric perturbations. 



6. CONCLUSIONS 

We described a numerical method based on a pseudo- 
spectral algorithm for studying the rapid variability properties 
of two-dimensional, viscous, hydrodynamic accretion disks. 
We demonstrated the ability of the spectral methods to han- 
dle correctly non-reflective boundary conditions and different 
stability conditions. We verified the implementation of the al- 
gorithm using various test problems. Also, we confirmed the 
non-axisymmetric instability of viscous spreading rings dis- 
covered by Speith and Kley. 

Spectral methods calculate the time evolution of a solution 
as an interaction of different waves. From a mathematical 
point of view, this approach naturally agrees with existing 
methods to study non-linear equations. It helps us to con- 
firm directly the result of linear mode analysis as well as tur- 
bulence theories. From a computational point of view, these 
methods are high order and produce accurate result in most 
problems. It is straightforward to extend the current algorithm 
to simulate 3D MHD accretion disks. In particular, because 
of the high order of the spectral methods, we can solve di- 
rectly for the vector potential and not for the magnetic field, 
in order to guarantee that the resulting field will be free of 
divergence at each time step. Also, spectral methods can be 
used to solve self-gravity problem without increasing signif- 
icantly the computation time. Finally, and most important, 
form a physical point of view, our algorithm provides another 
tool to set up numerical experiments and test astrophysical 
models. It is useful to confirm independently the results from 
other numerical simulations as well as test new models. 



APPENDIX 

A. CONVECTION-CONSERVATIVE MIXED FORMULISM 

As we describe in 33.21 we employ a convective-conservative mixed method. The non-linear term (v • V)v in the Euler equation 
is implemented in convective from, and all other terms are in conservation form. 



<9 r £: 



(Al) 
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a „ v d r (rT rr -rP) d^ P-t^ 
d t V r =-V r d r Vr- -^(d^Vr-Vj,) + — + -^ + — -^ 



o t v<t, = -v r d r v<f, {o^vj, + v r ) + — 



<> 






r, F _ d r {rEv r + rq r + rF r ) _ d^jEv^ + q^+F^) _ p . , , 



+ *!+< 
r r 



-2R. 



(A2) 
(A3) 
(A4) 



Here, the notation <9 r denote partial derivative with respect to time. The other notations d r and d^ denote the numerical derivative 
operators. Hence, the notation v r d r v r means, we first take a numerical derivative of the variable v r and multiply the result with v r , 
which is in convective form. On the other hand, d r (rEv r )/r means, we first use a temporary variable to store the produce rSv r , 
apply the procedure to calculate the numerical derivative on it, and then divide it by r. This is the conservative form. 
The viscosity tensor Ty in the above equation has the following general form 

(A5) 



Tij =2(fl r + fl s )eij+ I ^r + ^b-T/is I V-V. 

As we describe before, fi T , ^, and fi s are coefficient of radiative, bulk, and shearing viscosity. The strain rate tensor e,y written in 
cylindrical coordinate becomes 



e rr = d r v r 

1 ( r, Vd, dtpVr 
e-ri, =e<t,r=-z\ OrVj, + 

2 \ r r 

erhrh = 1 



and V • v in cylindrical coordinate is 

V-v = d,-v r + — + 
r 

For completeness, we provide again the viscous dissipation rate 



■ = 2(^ r + ^ s )(e, ) ) 2 + f /i r + /i b --A*s j (V-v) 2 



(A6) 

(A7) 

(A8) 
(A9) 

(A10) 



The functional forms of P, /i r , /i^, /U s , q, F, and F z can be change easily depend upon the physical assumption (where the defaults 
are described in ©. 
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